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Abstract. The Lieb- Robinson bound imphes that the unitary time evohition of an 
operator can be restricted to an effective hght cone for any Hamiltonian with short- 
range interactions. Here we present a very efficient renormahzation group algorithm 
based on this light cone structure to study the time evolution of prepared initial states 
in the thermodynamic limit in one-dimensional quantum systems. The algorithm does 
not require translational invariance and allows for an easy implementation of local 
conservation laws. We use the algorithm to investigate the relaxation dynamics of 
double occupancies in fermionic Hubbard models as well as a possible thermalization. 
For the integrable Hubbard model we find a pure power-law decay of the number of 
doubly occupied sites towards the value in the long-time limit while the decay becomes 
exponential when adding a nearest neighbor interaction. In accordance with the 
eigenstate thermalization hypothesis, the long-time limit is reasonably well described 
by a thermal average. We point out though that such a description naturally requires 
the use of negative temperatures. Finally, we study a doublon impurity in a Neel 
background and find that the excess charge and spin spread at different velocities, 
providing an example of spin-charge separation in a highly excited state. 
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1. Introduction 



Using ultracold atomic g quantum simulators, it has become possible to prepare 

states in almost perfectly isolated many-body systems and to monitor their time 
evolution [11 |2l [Sj SJ |5l |6]. At the same time, enormous progress in numerical 
renormahzation group methods has given us access to the dynamics of quantum models 
in one dimension (ID) [3 El El [ini EH US [13]. These algorithms are all based on 
approximating a quantum state as a matrix product in an optimally chosen truncated 
Hilbert space, an idea dating back to the density matrix renormahzation group (DMRG) 
by White [13]. This makes it now possible to study, both experimentally and numerically, 
fundamental questions about the relaxation dynamics and the role of conservation laws 
[H [15]. Furthermore, the applicability of the eigenstate thermalization hypothesis 
(ETH) — according to which each generic state of a closed quantum system already 
contains a thermal state which is revealed during unitary time evolution by dephasing 
[T6l [T71 [18] — can be investigated as well. 

In Sec. [2lwe present a new algorithm to study the unitary time evolution of an initial 
state in a ID quantum system. We concentrate on the case of a product initial state 
particularly relevant for experiment but note that the algorithm has been implemented 
also for thermal initial states. The main idea is to make use of the Lieb-Robinson 
bound [19] to efficiently simulate the system and to obtain results directly in the 
thermodynamic limit. Let us briefly recapitulate one of the main results of Refs. [T9l[20] 
which is the basis for our algorithm. We are interested in the time evolution of quantum 
systems starting from some initial state \^ i) where all connected correlation functions 
decay exponentially with a finite correlation length ^. The time evolution of a local 
operator acting on sites j,j + l,...,j + n can then be approximated by an 
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Figure 1. LCRG algorithm. Tlie Trotter-Suzuki decomposition of time evolution 
reveals the light cone structure. 



operator acting only in the effective light cone of the region [j, j + n] (see Fig. [T]) while 
being the identity operator outside of the light cone. More precisely, if o|^-^-^^](t) is the 
time evolved operator active only on sites which are at most distance / apart from the 
region [j, j + n] then 

\\%,j+n\it) - o\jj+n]it)\\ < const X cxp (^^- — (1) 

where vlr is the Lieb-Robinson velocity which is typically, in natural units, of 
the order of the interaction parameters of the model under consideration |20]. If 
vlrIA ^ I then the error of approximating oyj+„](t) by c>\j j+n]^'^) exponentially 
small. We show in Sec. [2] that a Trotter-Suzuki decomposition of unitary time evolution 
immediately leads to a light cone and that this light cone can be represented in a 
truncated Hilbert space using density matrix renormalization group (DMRG) techniques 
[llEIlESllMlEllESlEniiniETlEH]. in contrast to ground state and transfer matrix 
DMRG algorithms an explicit calculation of eigenvectors of the system, which is the 
computationally most costly step, is not necessary. This makes the new light cone 
renormalization group (LCRG) algorithm extremely fast and efficient. Furthermore, 
the implementation of local conservation laws — important for an effective numerical 
study — becomes particularly simple. The "speed of light" set by the Trotter- Suzuki 
decomposition is typically chosen to be much larger than the Lieb-Robinson speed 
vlr at which information spreads so that the algorithm directly yields results for the 
thermodynamic limit. Contrary to the infinite size time evolving block decimation 
(iTEBD) [13], however, it does not rely on translational invariance. In our paper we will 
demonstrate these advantages of the LCRG algorithm by studying several examples. At 
the same time we note that our approach does not solve the most fundamental problem 
of using matrix product states to investigate the time evolution of one dimensional 
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quantum systems: the linear growth of entanglement entropy with time, which restricts 
the applicability of such methods to the intermediate time dynamics. It can be shown 
under very general conditions that this is a fundamental property of unitary time 
evolution [191 ISO] which cannot easily be overcome. 

In Sec. [3] we will apply the LCRG algorithm to study the relaxation of a doublon 
lattice in ID fermionic Hubbard models. While the problem of a single doublon- 
holon pair has already been studied in ID [29J, our study is mainly motivated by 
the experimental and theoretical investigation of the decay of a macroscopic number 
of doublons in an ultracold fermionic gas on a three-dimensional optical lattice [5l |30] . 
First, we will present a test of the algorithm by studying the free fermion case where the 
time evolution can be calculated analytically. Next, we will investigate the differences 
in the relaxation dynamics between the interacting integrable and non-integrable cases 
as well as a possible thermalization in the long-time limit. In Sec. H] we will then 
demonstrate one of the major advantages of the LCRG algorithm: Even for systems 
without translational invariance, results in the thermodynamic limit can be obtained. 
In Sec. [5] we give a brief summary and an outlook on possible future applications of the 
algorithm. The supplementary material contains the executable code of the LCRG for 
the anisotropic Heisenberg model, which we have chosen as a simple example, as well 
as videos of the time evolution for the problem studied in Sec. HI 

2. The light cone renormalization group algorithm 

We present the LCRG algorithm to compute the time evolution 

(o[,,+„])^(t) ^ (M/,|e'^* 0[,,+„] e-'^*|vl/,) (2) 

of a local operator 0[jj+„] acting on sites j,j + + n. For the initial state 

= |si S2 . . .) we consider a product state with Sj denoting states in the local basis 
of dimension M. We note that with the help of ancilla sites also thermal states can be 
expressed using a product initial state followed by an imaginary time evolution, which 
is also performed using the light-cone algorithm. In this way we have implemented 
the real time evolution starting, e.g., from a highly entangled quantum state such as 
the ground state. We consider a Hamiltonian H = Y2j with nearest neighbor 

interaction; a Trotter-Suzuki decomposition of the unitary time evolution operator then 
leads to the 2D lattice shown graphically in Fig. [H It consists of local updates of 
two neighboring sites forward in time Tjj+i{6t) = exp{—ihjj+i6t) ("'|'" plaquettes) and 
backward in time Tjj^i{—6t) = rjj^-^^^St) ("|" plaquettes), where 6t is the Trotter-Suzuki 
time step. Unless there is an operator insertion, facing plaquettes trivialize and become 
the identity operator, 5t) rjj^i{6t) = 1 (shaded plaquettes). This yields the hght 

cone structure emanating from the local observable oyj^n] at time t. As long as the 
"speed of light" of the Trotter-Suzuki decomposition is larger than the Lieb-Robinson 
velocity vlr the expectation value is effectively evaluated in the thermodynamic 
limit. Neither translational invariance of the initial state nor of the Hamiltonian are 
required for this construction. 



CONTENTS 



5 




Figure 2. LCRG algorithm. Tlie liglit cone C grows with each time step by adding 
first a diagonal left transfer matrix L and then a diagonal right transfer matrix R. 

The LCRG algorithm is based on corner transfer matrices j3T| [32] to compute the 
growth of the hght cone with each successive time step 5t (Fig. [2]): the hght cone Ct at 
time t is multiphed from the left with the diagonal left transfer matrix L and then from 
the right with the diagonal right transfer matrix R to construct the new light cone for 
the next time step, Ct+st- Of course a direct implementation of this procedure would 
quickly break down because the Hilbert space of light cone states grows exponentially 
with time. Therefore, we use ideas from DMRG studies of dynamics in stochastic 
systems [25l [26] to represent both the light cone C and the transfer matrices L, R in 
a reduced Hilbert space of manageable dimension. A fully working implementation of 
this algorithm specialized to homogeneous systems in included in the supplementary 
material. 




Pl Plocal 



Figure 3. LCRG algorithm, (a)-(c) The light cone C grows to the left by contraction 
with a left transfer matrix L; the left and right transfer matrices L and R are augmented 
by a local plaquette T[5t). (d) The reduced density matrix is constructed from the 
forward and backward light cones by tracing only over the right site and block indices, 
(e) The local density matrix is obtained by tracing the forward and backward light 
cones over both left and right block indices. 
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In practice the time evolution proceeds in two half time steps (see Fig. In the 
first step, the light cone Cf[m;mr] has left and right block indices (Fig. |3]d) representing 
states in the (reduced) Hilbert space of dimension x, while the left transfer matrix 
Lt[miSimrSr] has again two block indices but also left and right site indices s;, Sr of 
dimension M (Fig. [3ti). Lt and Ct are contracted over their common block index to 
yield the new light cone Ct+st/2 half a time step ahead (Fig. Wp)'- 

Ct+5t/2[miSiSrmr] = ^ Lt[miSimSr] Ct[mmr] . (3) 

m 

The left transfer matrix Lt is enlarged by adding a plaquette at its upper right site index 
Sr (Fig.EJi): 

Lt+stMniiSis'irrirSrs'^] = ^ Lt[miSimrS] T[s'is'^sSr] . (4) 

s 

Similarly, a local plaquette is attached to the upper left corner of the right transfer 
matrix Rt to construct Rt+5t/2 (Fig. [3}:;). The initial conditions at t = are as follows: the 
block indices represent a single site with dimension mi = rrir = M, the intial light cone 
is Ct=o[mimr] = "^ilmirrir] for the product initial state |\I'/) on two neighboring sites, 
and the transfer matrices have the initial forms Lt=o[fniSimrSr] = ^/[^i'S] T[siSrSmr] 
and Rt=o[simiSrmr] = Y^s^l^t^rmis] "^ilsuir]. 

In order to bring C, L and R back into their original form the old block index m 
(dimension x) is combined with the adjacent site index s (dimension M) into a new 
block index m' = {ms) of dimension x' = Mx- The challenge is to limit the exponential 
growth of X with every time step. This is done by a renormalization step where a 
reduced density matrix is used to select the x most important basis states within the 
x'-dimensional Hilbert space. The reduced density matrix for the left block index is 
formed by combining the forward and backward light cones and tracing over the right 
site and block indices (Fig. [3li) 

pL[{ms'){ms)] = ^ Ct^st/2[{ms')srmr] Ct+5t/2[ims)srmr] (5) 

Srmr 

where we have used the fact that in unitary time evolution the backward light cone is 
the adjoint of the forward light cone. The reduced density matrix of dimension x' is 
by construction hermitean and has unit trace, pi is diagonalized, and the x states with 
the largest eigenvalues form the basis of the reduced Hilbert space. Optionally, one can 
retain all states such that the cumulative weight of the discarded states remains below a 
given threshold. We use a combination of both to obtain a reliable error control. Finally, 
the left block index of the light cone C and both block indices of the left transfer matrix 
L are projected onto this reduced basis, {misi) ^ mi. Analogously, the reduced density 
matrix p^ is formed by tracing over the left block indices to find a reduced basis for the 
right block indices, and subsequently the right block index of C and both block indices 
of R are projected onto the reduced Hilbert space. This completes the first half time 
step. 
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The second half of the algorithm works similarly by joining a right transfer matrix 
R to the right of the light cone (Fig. |2]), 

Ct+st[msiSrmr\ = ^ Ct+5t/2[mim] Rt+st/2[simsrmr] . (6) 

m 

At this stage the local density matrix 

piocaiit + 6t)[s'js'j^^SjSj+i]= ^ C*^st['i^is'js'j^^mr]Ct+st[miSjSj+imr] (7) 

minir 

is formed by contracting the forward and backward light cones over the left and right 
block indices, leaving open the site indices in the middle (Fig. [3^). The expectation 
value of a local operator oyj+n] is then obtained as 

{%,j+i]yit + 6t) = Tryj+i] (^piocai(t + St) 0[jj+i]) . (8) 

By multiplying further transfer matrices onto the left or right one can form also the 
local density matrix pyj+n] spanning more than two neighboring sites. For example, 
the density profile to the left of an impurity site is obtained by starting with pi and 
repeatedly multiplying L from the left onto the lower light cone and L* onto the upper 
light cone, until the desired distance from the impurity is reached. The remaining second 
half time step proceeds in complete analogy with the first part, growing L and R by one 
plaquette and renormalizing in turn the left and right block indices. 

Note that only summations and multiplications are required to build the light cone. 
This saves the most time-consuming step in standard transfer matrix DMRG algorithms, 
where one has to find the largest eigenvector of the transfer matrix. Only the density 
matrix p^ ^ has to be diagonalized, which dominates the computation time 0{M^x^)- 
Our algorithm therefore combines the speed of iTEBD [13] with the flexibility of TEED 
[12] to treat non-translationally invariant systems. Due to the local structure of the 
updates, conservation laws are easily implemented in our algorithm (see below). We 
note that instead of the first order Trotter-Suzuki decomposition shown in Fig. [1] also 
higher order decompositions can be easily implemented. 



3. Doublon decay in Hubbard models 

We use the LCRG algorithm with a second order Trotter-Suzuki decomposition to study 
dynamics in the ID fermionic Hubbard model 

Hu,v = - J (sVi+i.- + h.c.) + UJ2 (^n -\) Uji 

+ ^EK-1)K+i-1) (9) 

where J is the hopping amplitude, rij^^ = Cj ^Cj^„ and rij = rij^ + riji the occupation 
numbers, U the onsite, and V the nearest-neighbor potential. As initial states we will 
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consider the state where doubly occupied and empty sites alternate, and the Neel 

state, \^n)- These state are given explicitly by 

i^d) = n ' 1^^) = n 4-+it4-jo) (10) 

3 j 

where |0) denotes the vacuum. For these states, we want to investigate the time 
dependence of double occupancies, dj = nji^nji, the staggered magnetization, mj = 
{—lySj = {—ly^rij-f — nji)/2, and of the operator Wj = (— 1)%^, measuring the charge 
imbalance between even and odd sites. Before discussing the numerical results, we first 
want to establish a number of relations between these three operators in the case where 
the nearest neighbor repulsion vanishes, = 0. The model with V ^ will be studied 
in Sec. KM 



3.1. Duality relations for the integrahle model 

For V = Q the model ([9]) becomes the integrable Hubbard model. Apart from the 
special symmetries responsible for the integrability of the model by Bethe ansatz, there 
are other symmetries in this case which allow us to establish various relations between 
the states and operators: 
(a) There is a unitary duality transformation 

w=n(^^t+(-i)sv) (11) 

3 

relating the repulsive [U > 0) and the attractive {U < 0) Hubbard models. This 
transformation leads to Wcj^^U = {—lYcj^, Wcj\JA = Cji so that the kinetic energy part 
in Eq. stays invariant while U — )■ —U. For the operators we find 

dj ^nji-dj, (12) 
m,^ (-1)^(1 -n,)/2, (13) 

^ (-l)'(l + ^.i-^.t)- (14) 
For the initial state it follows that U^'^d) = j\f) , U^"^ j\f) = I^^d), assuming an even 
number of lattice sites L. For the expectation values (see also Eq. ([2])) 

L 

oh{t)^J2(o,)Ut)/L, (15) 

the duality transformation implies 

d^{t) = 1/2 -d%{t), (16) 

^^it)= -nJ%{t)/2, (17) 

^uit)= -w%{t)/2^Q. (18) 

The expectation values in the last equation (|T8|) have to vanish identically because of 
the particle-hole (spin inversion) symmetry of the initial states, respectively. The second 
identity (fT7|) . furthermore, shows that the decay of the staggered magnetization can be 
studied in a realization of a fermionic Hubbard in cold atomic gases without the need 



CONTENTS 



9 




Figure 4. (a) it) for free SFF (circles) and for free SLF (squares) — note tliat tlie 
time scale in [33] is stretched by a factor 2. (b) d^{t) for free SFF with x as indicated. 
In both cases 5t = 0.1 and lines denote the exact results, (c) Absolute error 5d^ [t) for 
free SFF with x 5000 and 5t = 0.2,0.1,0.05,0.02 (in arrow direction), (d) S'c„t(t) 
for free SFF with 6t — and x a-s indicated. 



to address the spin degree of freedom directly in a measurement. 

(b) On a bipartite lattice, A ® we can furthermore apply the transformation 
Cja — 7- icjcr for i E A (j G -B), respectively. This leads to J — t- —J, U ^ U and 
therefore H_jj — )■ —Hu. This results in d^{t) = 1/2 — d^{—t) and similarly for the 
other identities. 

(c) Finally, we can use the time reversal invariance of the expectation values. Using all 
three symmetries we find 

dEit) =d%{t) = 1/2 -d^it), (19) 
m^{t) = m%{t) = -wE{t)/2. (20) 

The relaxation dynamics we will consider here is therefore independent of the sign of U 
and the same information is obtained by starting either from \^ d) oi \^ n)- 

3.2. Testing the LCRG algorithm: The free fermion case 

To test the LCRG algorithm, we first study the free spinful fermion (SFF) case t/ = 
where the dynamics can be calculated exactly. We find d^^Q(t) = (1 + Jo(4Jt))/4 and 
''T^u=oi^) ~ -^o(4Jt)/2 with Jo the Bessel function of the first kind. In the free SFF 
case, the dynamics of electrons with spin up and spin down is completely decoupled. 
Therefore, we can also use free spinless fermions (SLF) to calculate ttlq {t) with a spinless 
particle representing either the presence of a spin up or a spin down. Then we need to 
keep only ^/x states to simulate the dynamics with the same accuracy. In Fig. HJ^a) the 
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LCRG results for itiq (t) for free SFF and SLF are compared to the exact result. For 
free SLF with x = 20000 block states we are able resolve 6.5 oscillations compared to 
the 5 oscillations which have been resolved in |33] by iTEBD. We emphasize that for the 
Hubbard model (M = 4) with the conservation laws for spin and charge implemented 
and X = 2000 states kept, each time step takes only ~ 30 seconds on a standard PC 
without parallelization. This is 260 x faster and uses 12 x less memory than without 
conservation laws, because the largest diagonal block of the reduced density matrix 
Pl,r has only 200 states. For SLF (M = 2) the speedup is still 40 x with 5x less 
memory and a largest block of 450 states. In Fig. Hl^b) the results for are shown 

where x is varied. The error of the simulation up to tmax where the simulation starts to 
deviate from the exact result is completely dominated by the error of the Trotter-Suzuki 
decomposition (see Fig. ID^c)) and is of order for the second order decomposition 
used here. Importantly, tmax is determined only by x ^md results with in principle 
arbitrary accuracy can be obtained for t G [0, tmax] by reducing the time step 6t or 
using a higher order Trotter- Suzuki decomposition, since the number of RG steps is not 
restricted. 

The algorithm breaks down when the spectrum of the reduced density matrix 
Ps = Pl,r becomes dense. A suitable measure is the entanglement entropy 

^ent(t) = -Trp.lnp, < In(Mx) (21) 

with Mx = dimp^. The entanglement entropy is shown in Fig. 11(d) and increases 
linearly with time. We want to remind the reader once more that the linear increase of 
the entanglement entropy seems to be a fundamental property of unitary time evolution 
[20] which cannot easily be overcome and limits the simulation time. 

The LCRG algorithm actually does provide an intuitive picture for this behavior: 
Facing plaquettes outside of the light cone (shown shaded in Fig. [T]) trivialize, thereby 
connecting a local degree of freedom at the edge of the lower light cone with one on 
the upper light cone by a Kronecker delta. The number of these Kronecker delta bonds 
between the lower and upper light cone increases linearly with time and determines the 
entanglement entropy between the light cones. This is very similar to the entanglement 
entropy of a spin-1/2 Heisenberg chain: In this case the ground state can be represented 
in a resonating valence bond (RVB) basis. If the chain is now split into two semi- 
infinite segments then the entanglement entropy has been shown to be proportional 
to the number of RVB bonds connecting the segments [3l]- A breakdown of the 
simulation is observable as a deviation from the linear growth of Sent and occurs when the 
entanglement entropy is close to the bound, Sentit) ~ In(Mx), i.e., when all eigenvalues 
of ps have comparable magnitude thus making further RG steps impossible. 

3.3. Results for the Hubbard model 

Next, we study d^{t) in the interacting Hubbard model, a situation which can be realized 
in ultracold gases [5]. For times Jt ^ min{J/|?7|, 1} the relaxation is independent of 
the interaction strength and follows the short-time expansion of the free fermion result 
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lUI/J lUI/J 

Figure 5. Evolution of double occupancy in Hubbard model, (a) d^{t) with 
X — 20000, 5t — 0.1 (symbols), and fits (lines), see text, (b) Slope of the entanglement 
entropy, 5cnt ^ aJt. (c) Fitted exponent a of the power-law decay, see Eq. p2|) . 

(i^(t) ~ 1/2 — 2(Jt)^, see Fig. [5]^a). Thus, in order to see the effect of interactions, 
systems at times Jt ^ have to be studied. In units of the hopping amplitude 

J we can simulate longer times the larger U is. This is a consequence of the slower 
increase of Senti^) ~ o-Jt as shown in Fig. [5t^b). For large U we find that the slope of the 
entanglement entropy is given by a ~ J /\U\., i.e., the simulation time is proportional to 
^max ~ \U\/ J"^ and therefore set by the inverse of the magnetic superexchange interaction 
~ J^/|f/|. It is clear that the slope of the entanglement growth becomes smaller the 
closer the initial state is to an eigenstate of the Hamiltonian: the eigenstate stays 
invariant under time evolution and no additional entanglement entropy is generated. 
Comparatively long times can therefore be simulated, in particular, if the time evolution 
of the ground state with a weak perturbation is studied as, for example, in Ref. [29]. 

At times Jt ^ 1 the relaxation in the free SFF case is given by = 
(1 + Jl{AJt))/A ~ [1 + (47rt)-i(l + cos(8Jt - 7r/2))]/4. This motivates us to fit the 
time dependence at finite U by the function 

d^if) = rf^(oo) + e-""[A + 13 cos{Qt - 0)]/t" (22) 

in the regime 1.5 < Jt < Jt^s.x- Such fits are shown as solid lines in Fig. [5]^a_). In 
all cases 7 < 10^^, i.e., we do not find evidence for a finite relaxation rate 7 A 
relaxation following a power law has also been observed at intermediate times in a ID 
Bose-Hubbard model starting from an initial state with one boson on every second site 
[35]. On the other hand, the relaxation for the XX Z model when starting from a Neel 
state has been interpreted in terms of an exponential decay [33]. Our fits point to a pure 
power-law decay with an exponent a which increases with increasing \U\, see Fig. [5]^c). 
The asymptotic value d^{oo) also increases and reaches 1/2 in the limit \U\ — > 00, see 

I We note that the fits for large \U\ are more ambiguous because d^{oo) is reached more quickly. 
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the y = data in Fig. [Tl^b) below. We emphasize that d^{t) = d^jj(t), i.e., repulsive 
interactions lead to a binding of doublons the same way as attractive interactions do |36] . 
For doublons moving on a 3D lattice it has been argued that for repulsive interactions, 
U > 0, much larger than the bandwidth, many-body scattering processes are needed 
to dissipate the doublon energy, which leads to an exponentially small relaxation rate 
7 ~ exp(— f/ / J) [30], |5]. In our simulations we do not see indications for a corresponding 
crossover time scale ~ I/7 at which exponential relaxation might set in. The power law 
decay in the Hubbard model (or, at least, the very small relaxation rate) might be 
a consequence of the infinitely many local conservation laws leading to integrability. 
It is important to stress, however, that our numerical data for the intermediate time 
dynamics cannot finally resolve the question whether or not exponential relaxation 
does exist. In the next paragraph we will, however, give further support that the 
fits with Eq. ( l22l) describe the relaxation at long times correctly by showing that the 
the asymptotic value d^{oo) obtained from the fits agrees very well with a thermal 
expectation value. 



3.3.1. Long-time limit and thermalization According to the eigenstate thermalization 
hypothesis (ETH) [161 each initial state — which can be represented as a 

superposition of eigenstates of the Hamiltonian — already contains a thermal state. This 
thermal state is revealed during time evolution due to dephasing effects between the 
different eigenstates. We say that the system has thermalized if the long-time average 

lim - I dt{^ i\o{t)\<^ i) (23) 



r— >oo T 

is equal to the thermal average {o)\^ in an appropriately chosen ensemble with the 
intensive variables Aj. Note that this definition only demands that o = {o)\^, i.e., time 
dependent fluctuations in (\l/7|o(t)|\&/) can, in principle, remain large even for t — )• 00. 
This will, in particular, be true for free models where no relaxation mechanisms exist and 
the concept of thermalization therefore has limited meaning. In interacting models, on 
the other hand, we expect that o{t — )■ 00) = o, i.e., time-dependent fluctuations vanish 
in the long-time limit. In this case, the system is expected to have truly thermalized if 
o(t -> 00) = {o)x,. 

The appropriate thermal ensemble is determined by the set of conserved quantities 
Qi with [H,Qi] = 0. Obviously, {'^i\Qi(t)\'^j) = const which means that the intensive 
variables (Lagrange multipliers) Aj have to be determined such that 

(^/IQil^/) = {Q^)xr (24) 

Every quantum system has two types of conserved quantities: local and non-local. We 
call a conserved quantity local if it can be expressed as a sum of local densities acting 
only on a finite number of lattice sites|§| For a generic system, usually only very few 
local conserved quantities such as the particle number operator and the Hamiltonian 



§ Equivalently, a conserved quantity in a field theory is local if it can be written as an integral of a 
fully local operator density. 
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itself exist. Any quantum system in the thermodynamic hmit has, on the other hand, 
infinitely many non-local conserved quantities, as, for example, the projection operators 
\En){En\ onto the eigenstates of the system. It is not clear if these non-local conserved 
quantities play any role in determining thermalization or transport [371 [TOt [TT] . In 
studies of thermalization they are usually simply neglected. 

The ID Hubbard model is integrable by Bethe ansatz and has infinitely many local 
conserved quantities which can be constructed explicitly from a family of commuting 
transfer matrices [3H1I39]. In this case, one should, in principle, consider a generalized 
Gibbs ensemble with a Lagrange multiplier Aj for each local conservation law. However, 
such ensembles are impossible to handle in the thermodynamic limit except for the 
simplest free particle models [18]. Here we will instead consider the usual canonical 
ensemble. The effective temperature Teg then acts as a Lagrange multiplier determined 
such that 

{^D\H\^D)/L = ^Ti{He-'"^'^^] (25) 

where Z = Tre~^/^'=* is the partition function and the energy = U/A 

is conserved during time evolution. Since the spectrum of eigenenergies per site is 
bounded for a lattice model, the use of negative temperatures is natural with 1/T — )■ 0^ 
corresponding to the case of maximum entropy. In the following we denote the thermal 
average in the canonical ensemble by {d)u^T and it is easy to see that {d)u,T = {d)-u-T 
holds. The duality transformation ( ITTl) . ( IT6l) furthermore implies {d)u,T + {d)-u,T = 1/2 
for all T. In particular, ((i)f/>o,T>o < 1/4 and ((i)(/<o,T>o > 1/4 with both being equal to 
1/4 in the limit T — )■ oo. We calculate the thermal average {d)u,T using a transfer-matrix 
DMRG algorithm [211 1221 SHI HI] and find that Teg is negative (positive) for the repulsive 
(attractive) model, respectively. The dependence of the effective temperature Tes on 
interaction strength is shown in Fig. Wi^) together with the results for the extended 
Hubbard model which are discussed in the next subsection. The comparison between 
the double occupancy extrapolated in time, d§{oc), and the thermal double occupancy 
{d)T^f^ shown in Fig. Wljo) yields excellent agreement. This seems to suggest, on the 
one hand, that the extrapolation using Eq. (l22l) is appropriate leaving very little room 
for an additional exponential decay which possibly could set in at a longer time scale. 
On the other hand, it also seems to suggest that the other local conservation laws of 
the Hubbard model have very little influence on the relaxation of double occupancies. 
To qualitatively understand the latter property we can think of writing the operator 
d as a sum of projections onto all the conserved quantities which form a basis of 
the operator space [371 [IH]- At least for large U it is then clear that the projection 
onto the Hamiltonian, which does contain the operator d itself, will give the dominant 
contribution to the thermal expectation value. 

3.4- Results for the extended Hubbard model 

Finally, we consider the non-integrable extended Hubbard model obtained by turning 
on the next-nearest neighbor repulsion V in Eq. ([9]). The duality relations used for the 
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Hubbard model are then no longer valid because V ^j{nj — l){nj^i — l) — )■ AV S^^S^^i 
under the duality tranformation, Eq. ffTTj) . However, we still have d^y{t) = d^jj_y{t) 
for the time evolution and {d)u,v,T = (d) -u,-v,-t for the thermal expectation value since 
these relations only rely on the lattice being bipartite and time reversal symmetry. We 
focus in the following on V < U/2 corresponding to a spin-density wave state in the 
ground state phase diagram |10l HI]. We note that for V > U/2, |\I/d) is close to the 
charge- density wave (CDW) ground state and long simulations in time are possible with 
d-uvi^) staying close to 1/2 and showing revival oscillations (data not shown). For the 
case V = U/ 2 — which is approximately at the phase transition line from the spin-density 
to a charge-density wave state [iQl SI] — we find a qualitatively different behavior than 
in the Hubbard model (see Fig. [6]). Here d^y{t) decreases at long times with increasing 



0.4r-pi 1 1 1 1 1 1 1 1 1 1 1 1 1 1 — I — r 




Figure 6. Extended Hubbard model: d^y{t) for V = U/2 with x = 10000, St = 0.1 
(symbols), and fits (lines). Inset: Relaxation rate 7 extracted from the fits. 

interaction strength. In general, d^y{t) can increase or decrease at long times with 
increasing interaction strength depending on the ratio V/U . Using the same fit function 
( !22l) as before we find that it is no longer possible to describe the relaxation dynamics 
by a pure power law decay. Instead, we now find a finite relaxation rate 7 as shown in 
the inset of Fig. |6l 

3.4- i- Long-time limit and thermalization Investigating again a possible thermalization 
we can shed some light on the observed dependence of the extrapolated value d^y{oo) 
on the ratio V/U . For V ^ the model is no longer integrable and — if thermalization 
does occur — the final state should be fully described by the canonical ensemble. The 
energy during the time evolution is now fixed to 

{^d\H\^d)/L = ^-V (26) 

and determines via the relation fl25]) the effective temperature shown in Fig. W^a). For 
V = U/4: it follows that l/Tgfj = and {d)ify=u/4 = 1/4. In the repulsive case, U > 0, 
the energy is negative for V > U/A leading to Tgg > 0. For V < f//4, on the 
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U/J lUl/J 



Figure 7. (a) Inverse effective temperatures for different ratios of V/U with [/ > 0. 
Toff changes sign for C/ < 0. (b) d^ y{(X)) (filled symbols connected by solid lines) 
compared to the thermal average ((i)T„tf (open symbols connected by dashed lines). 



other hand, the energy is positive and therefore Teg < 0. For f/ < the signs of Tgfj 
are reversed. Using again a transfer matrix DMRG algorithm to calculate the thermal 
expectation value {d)^^^ at temperatures Teg we can compare with the extrapolated 
value (i^y(oo) from the time evolution, see Fig. [D^b). Compared to the pure Hubbard 
model the agreement is not quite as good and the deviations increase the closer the 
ratio of the interactions is to the critical line V = U /2 and also the larger U is. This 
does suggest — assuming that the system will finally thermalize — that the numerically 
obtained intermediate time dynamics is not sufficient to fully extract the long time 
behavior. A possible explanation is that two different relaxation processes exist in this 
case: a fast one at short time scales leading to a pre-thermalized state and a slower 
one setting in at Jt ^ exp(f//J) [301 IS]- This might also explain the non-monotonic 
behavior of 7(f/, V/U = const) obtained when fitting with a single relaxation rate, see 
inset of Fig. [61 

4. Application to a non-translationally invariant case 

On of the main advantages of the LCRG algorithm is that it allows to study the time 
evolution of one-dimensional quantum systems in the thermodynamic limit even if the 
initial state and/or the Hamiltonian are not translationally invariant. As an example, we 
consider the time evolution in the = Hubbard model ([9|) of the non-translationally 
invariant state 

1^^) = cj^l^^v) = 4^ n 4,+it4,j0) (27) 

j 

obtained by adding an additional electron at site j = to the Neel state, Eq. (11 01) . In 
Fig. [HI LCRG results for the dynamics of the excess charge and the excess spin defined 
by 

(^r) = - (^?) ; (^D = (^.) - (^?) (28) 

are shown for U/J = 2. Here (n^^) = 1 ((s^^)) are the background charge (spin) 
densities, respectively, obtained from the time evolution of {"^n), i-e., from a system 
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Figure 8. Time evolution starting from the initial state IvPat), Eq. p7| . for U/ J = 2 
with X = 1024 states kept. Shown are results for times Jt = 3, 2.4, 1.8, • • • ,0 (from 
bottom to top) for (a) the excess charge {nj^^), and (b) the excess spin (s^™). 
Subsequent curves are shifted by 0.2 for clarity of presentation. The dashed lines 
connect the points Xc(s) , see Fig. HI 

without the additional electron. Videos of the time evolution for several other interaction 
strengths U/J are presented in the supplementary material. We observe that (nj^^) and 
{sY'^) spread out into the lattice with different velocities clearly revealing the light- 
cone structure. For the non-translationally invariant problem considered here we have 
not implemented the conservation laws yet and the results shown in Fig. [8] have been 
obtained by keeping a relatively moderate number of states, x = 1024. Although the 
simulation time is therefore smaller than in the translationally invariant cases discussed 
in the previous chapters, we want to stress that the evolution of the Neel background 
is simulated in the thermodynamic limit and thus very different from that in a small 
system tractable, for example, by exact diagonalization. 

Different charge and spin velocities for the one-dimensional Hubbard model starting 
from an initial non-equilibrium state have already been observed in Ref. [12]. In this 
case, an initial state was considered where the ground state was perturbed by a small 
local charge and spin imbalance thus allowing to extract the velocities of the elementary 
spin and charge excitations which can also be obtained by Bethe ansatz. Our initial 
state, on the other hand, is a highly excited state and the charge and spin velocities are 
not related to those of the elementary excitations. 

To extract the charge velocity Vc and the spin velocity Vg for our initial state we 
take the point Xc{s)it) where the tail has reached half of the height of the first peak of 
the charge (spin) distribution as reference point. Fig. [9](a) shows that Xc(s)i't) depends 
linearly on time. The dependence of the velocities on the interaction strength U/J is 
shown in Fig. [9](b). We find a charge velocity Vc ~ 2 J independent of U. For U = the 
charge and spin distributions are identical for all times and Vg = Vc, but for increasing U 
the spin velocity Vs decreases. These results can be qualitatively understood as follows: 
the excess charge sees a uniformly charged background and therefore moves unimpeded 
with the Fermi velocity of the non-interacting system, Vc ~ vp = 2 J. The excess 
magnetization, on the other hand, moves with a velocity proportional to the effective 
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Figure 9. (a) Time evolution of the reference points Xc[s) ^.s defined in the text for 
the excess charge (circles) and for the excess spin (squares) with U/J ~ 2. The lines 
are linear fits, (b) Charge and spin velocities as a function oiU/ J. 

spin superexchange ~ J^/f/ which therefore decreases with increasing U . For large U 
we find, furthermore, that the spin dynamics becomes more comphcated. While most 
of the excess magnetization remains inert in this limit, an additional small staggered 
part appears which spreads out with a velocity close to the charge velocity (data not 
shown) . 

5. Conclusions 

The development of new spectroscopic techniques to study ultracold quantum gases on 
optical lattices with very good spatial and time resolution has put the topic of non- 
equilibrium dynamics in quantum systems firmly back onto the agenda. Particularly 
interesting from a fundamental persepective is the dynamics in one dimension where 
many of the standard lattice models such as the Heisenberg, the t — J, and the Hubbard 
model are integrable, i.e., these systems have an infinite number of local conservation 
laws. Since a conserved quantity stays invariant under time evolution, the presence of 
many conserved quantities is expected to severely restrict the dynamics of the quantum 
system as a whole and might even prevent the system from reaching thermal equilibrium. 
To investigate such questions numerically, different algorithms have been developed: The 
time-dependent density-matrix renormalization group (t-DMRG) and the time-evolving 
block decimation (TEBD), in particular, allow one to access the intermediate time 
dynamics by representing the time evolved state as a matrix product state. Among 
the aims when developing numerical algorithms to study the time evolution are (a) 
the thermodynamic limit, (b) the flexibility to simulate different systems, (c) a high 
computational efficiency, and (d) to simulate the system for as long as possible. 

Here we have presented a new algorithm which does make progress concerning 
many of the points mentioned above. The light cone renormalization group (LCRG) 
algorithm is highly efficient by simulating at each time step only that part of the lattice 
(a light cone) which is affected by the time evolution. The results obtained are directly 
for the thermodynamic limit. Contrary to the infinite size TEBD, the LCRG algorithm 
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does not rely on translational invariance. It is therefore extremely flexible and can, in 
particular, also deal with non-translationally invariant problems. It cannot, however, 
simulate the quantum system for times significantly longer than other algorithms based 
on matrix product states. We have shown that the lattice path integral representation of 
unitary time evolution — forming the basis for the LCRG algorithm — provides a simple 
picture for the linear entanglement growth with time which restricts the simulation time 
of such algorithms. An executable sample code for the anisotropic Heisenberg model is 
provided in the supplementary material. 

We used the LCRG algorithm to study quench dynamics in the integrable Hubbard 
model starting from a state with every second site doubly occupied and found indications 
for a pure power-law relaxation. For the extended non-integrable Hubbard model, on 
the contrary, the relaxation appears to be exponential. In both cases we found that the 
time evolved state in the long-time limit seems to be close to a thermal state supporting 
the eigenstate thermalization hypothesis. 

Finally, we demonstrated that the LCRG algorithm can also be used to study the 
time evolution of non-translationally invariant initial states. For the Neel state with 
one site occupied by a doublon we showed that the excess charge and excess spin spread 
with finite, but different, velocities. We extracted the dependence of the velocities on 
the strength of the Hubbard interaction U which can be used for comparison with future 
experiments. Videos of the time evolution for different parameter sets can be found in 
the supplementary material. 

As an outlook we want to emphasize that future applications of the LCRG algorithm 
to other non-translationally invariant setups such as impurity problems, disorder, and 
to systems with trapping potentials are feasible. 
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